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Abstract. In this paper, we design high order accurate and stable finite dif- 
ference schemes for the initial-boundary value problem, associated with the 
magnetic induction equation with resistivity. We use Summation-By-Parts 
^ (SBP) finite difference operators to approximate spatial derivatives and a Si- 

multaneous Approximation Term (SAT) technique for implementing boundary 
conditions. The resulting schemes are shown to be energy stable. Various nu- 
merical experiments demonstrating both the stability and the high order of 
accuracy of the schemes are presented. 

< 

^ 1. Introduction 

^ Many interesting problems in astrophysics and engineering involve evolution 

y of macroscopic plasmas, modeled by the equations of MagnetoHydroDynamics 

(MHD). These equations ([ ]) are a system of convection-diffusion equations with 
the magnetic resistivity and heat conduction playing the role of diffusion. Many 
^ applications like plasma thrusters for deep space propulsion and electromagnetic 

pulse devices ([ ;]) involve small (but non-zero) values of the magnetic resistivity. 
Hence, the design of efficient numerical methods for the resistive MHD equations 
is essential for simulating some of the afore mentioned models. 

Numerical study of the ideal MHD equations (where magnetic resistivity and 
other diffusions terms are neglected) has witnessed considerable progress in recent 
I years and a variety of numerical methods are available (see [(i] for a review of 

the available literature). The design of numerical schemes for the resistive MHD 
equations has not reached the same stage of maturity as the presence of magnetic 
. ^ resistivity complicates the design of stable methods even further. Given the formi- 

dable difhculties, study of prototypical sub-models (that mirror some, but not all of 
^ the difficulties of the resistive MHD equations) can be a useful guide for obtaining 

robust methods for the resistive MHD equations. 

In this paper, we consider the magnetic induction equations with resistivity. Re- 
cent papers ([2, ')]) have pointed out the role that the magnetic induction equation 
(without resistivity) plays in the design of numerical schemes for the ideal (inviscid) 
MHD equations. The induction equations with resistivity can play a similar role 
for designing stable methods for the resistive MHD equations. Our goal in this pa- 
per is to design stable and high-order accurate numerical schemes for the magnetic 
induction equations with resistivity. 
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We start with a brief description of how the equations are derived. In a moving 
medium, the time rate of change of the magnetic flux across a given surface S 
bounded by curve dS is given by (see [ i ] ) : 

^ jB-dS = J^-dS + j)Bxu-dl + J (div(B))u ■ dS + e ^ 3 ■ dl, 

s s as s OS 

where the unknown B = B(x, G M'^ denotes the magnetic field, J = J(x, e M.^ 
the current density and x = {x,y,z) are the spatial coordinates. The current 
density is given by: J = curl(B). The parameter e denotes the magnetic resistivity, 
and u(x, t) the (given) velocity field. 
Using Faraday's law: 

(1.1) -^y"B-dS = ^E'-d/, 

s as 

Stokes' theorem, the fact that the electric field E' = in a co-moving frame and 
E' = E + u X B we obtain, 

(1.2) — + curl(B X u) = -udiv(B) - ecurl(curl(B)). 
at 

Magnetic monopoles have never been observed in nature. As a consequence, the 
magnetic field is always assumed to be divergence free, i.e., div(B) — 0. Using this 
constraing in (1.2), we obtain the system: 

dtB + curl(B X u) = -ecurl(curl(B)), 

^^•^^ div(B) = 0. 

The above equation is an example of a convection-diffusion equation. The version 
obtained by taking zero resistivity (e = 0) in (1.3) is termed the magnetic induction 
equation ([II]). A standard way to obtain a bound on the solutions of convection- 
diffusion equations like (1.3) is to use the energy method. However (1.3) is not 
symmetrizable. Consequently it may not be possible to obtain an energy estimate 
for this system. 

On the other hand, (1.2) is symmetrizable. We use the following vector identity 
curl(B X u) = Bdivu - udiv(B) + (u • V) B - (B • V) u 

= (u'B)^ + (u^B)^ + (u'B)^ ~ udiv(B) - (B • V)u, 

and rewrite (1.2) in the form, 

dtB -I- (u • V) B = -B(divu) + (B • V)u - ecurl(curl(B)) 

^^'"^^ =M(i:>u)B-ecurl(curl(B)), 

where the matrix M{Du) is given by 

?yU^ — dzU^ dyU^ dzU^ 

M{Du) = I d^u^ -d^u^ - dzv? dzU^ 



dxU^ dyV? —dxU^ — dyU^ 



Introducing the matrix, 



C = — I dxu'^ dyU^ dzU^ 



JyL 



dzu\ 
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(1.2) can also be written in the following form, 

(1.5) dtB + 9^ (A^B) + dy (A^B) + (A^B) + CB = -ecurl(curi(B)), 

where A' — u^I for i — 1,2,3. Note that the symmetrized matrices in (1.5) are 
diagonal and that the coupling in the equations are through both the lower order 
source terms and the viscous terms. 

Furthermore, taking the divergence of both sides of (1.2) we obtain, 

(1.6) (div(B))t + div (udiv(B)) = 0. 

Hence, if div(B(x,0)) = 0, it follows that div(B(x,t)) = for t > 0. This implies 
that all the above forms (1.5), (1.3) and (1.2) are equivalent (at least for smooth 
solutions) . 

Although the magnetic induction equations with resistivity are linear, the coef- 
ficients are functions of x and t. Therefore, closed form solutions are not available 
and we must resort to numerical methods in order to calculate (approximate) solu- 
tions. Consequently, it is important to design efhcient numerical methods for these 
equations. 

As mentioned before the magnetic induction equations is a sub-model in the 
resistive MHD equations. Hence, design of stable and high-order accurate numerical 
schemes for the viscous induction equations can lead to robust schemes for the non- 
linear resistive MHD equations. 

The presence of the div(B) = constraint leads to numerical difficulties. Small 
divergence errors may change the nature of results from numerical simulations (see 
[4, 12] for details on the role of divergence in ideal MHD codes). Our approach 
to treating the constraint follows the method developed in [7, ] and involves 
discretizing (1.5). TFurthermore, a proper discretization of the symmetric form 
(1.5) yields energy estimates. These estimates are vital in proving existence of weak 
solutions. We will approximate spatial derivatives by second and fourth order SEP 
("Summation By Parts") operators. The boundary conditions of both the Dirichlet 
and mixed type are weakly imposed by using a SAT ( "Simultaneous Approximation 
Term"). This work is an extension of the SBP-SAT schemes for the case without 
resistivity (e — 0) found in a recent paper [ ]. 

We would like to emphasize that other numerical frameworks like mixed finite 
elements, discrete duality finite volume or mimetic finite differences might also 
lead to stable schemes for approximating the induction equations with resistivity. 
However, we are not aware of any papers that have approximated the resistive 
induction equations with these approaches. 

The rest of this paper is organized as follows: In Section 2, we state the energy 
estimate for the initial-boundary value problem corresponding to (1.4) in order to 
motivate the proof of stability for the scheme. This is done for mixed type and 
Dirichlet boundary conditions. In Section 3, we present the SBP-SAT scheme and 
show its stability with both Dirichlet and mixed boundary conditions. Numerical 
experiments are presented in Section 5 and conclusions from this paper are drawn 
in Section 6. 

2. The Continuous problem 

For simplicity and notational convenience, we restrict ourselves to two spatial di- 
mensions in the remainder of this paper. Extending our results to three dimensions 
is straightforward. 
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In two dimensions, (1.4) reads 

Bj + AiB^ + AsBj^ 

(2.1) 



where - V x (V x B) = 



CB = -eV X (V X B) , 



xy) 



and 



Ai = 



A2 = 



w u 
w 

T 



c = 



dyU^ 

— drU 



(2.4) 



curl(B) X n 



and 



with B = {B\B^) and u — denoting the magnetic and velocity fields 

respectively. Throughout this paper, we consider (2.1) in a smooth domain Q. 
One can extend our results to general piecewise smooth boundaries by a standard 
procedure. We augment (2.1) with initial conditions, 

(2.2) B(x,0) = Bo(x), xe^!, 

and Dirichlet or mixed boundary conditions with homogeneous boundary data. The 
Dirichlet boundary conditions are given as 

(2.3) B(x,i)=0, xedfl. 

In order to specify the mixed boundary conditions, we need some notation. Let 
n(x) denote the outward pointing unit normal at a point x G dQ. Define 

-n^{Bl-Bir 
^ n\Bl-Bl) 

(2.5) ^\iB) = Bl-Bl 

for B = {B^,B^) and n = {71-^,71^). Furthermore, let Qriin denote the part of d^l 
where the characteristics are incoming, i.e., 

d^in = {x e 917 I n(x) • u(x) < 0} . 

With this notation, the mixed boundary conditions read 

/3(n(x) • u(x)) B(x, t) + c^l(B(x, t)) x n(x) = 0, for x G dilin, 

cml(B(x,i)) X n(x) = 0, for x e dn\ dfli^, 

where /? < — l/(2e) is a given number. 

In order to motivate the complicated calculations required to show stability in 
the discrete case, we start by explaining how stability is proved in the continuous 
case. We assume that the solution, B, is sufficiently regular for our calculations to 
make sense. 

Theorem 2.1. Lei B(x, <:) be a solution of the problem (2.1) and (2.2) with bound- 
ary conditions (2.3) or (2.6). There exist positive constants a (depending on u and 
its first derivatives) and K , such that 

2 



(2.6) 



m;t)\\ 



L2(o) 



L2(o) 



< Ke° 



!Bo| 



L2(o) • 



curl(B(-,t)) 

Proof. Multiplying (2.1) by B^ and then integrating in space, we get 

/ B'^dtBdxdy + [ (B'^AiB^ + B'^AzB^ - B'^CB) dxdy 
Jn Jn 
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-e / (V X (V X B)) dxdy, 



which imphes 
2dt J a 



2 

so that, 

ld_ 
2dt 



dxdy + ^ y dxdy 
^ / B'^ {2C + div(u)) B dxdy ~-[ {u ■ n) ds + e [ (b • {c^\{B) x n) 



B2 

n 



< a 



dxdy ^" ^ y ^curl(-^)) dxdy 

I B^ dxdy - \ ( B^ (u • n) rfs + e / (b • (c^l(B) x n)) ds. 
Jn 2 Jq^ Jqq \ ) 



From the above relation, we see that applying Dirichlet boundary conditions, (2.3), 
and integrating in time gives the required result. For the mixed boundary conditions 
(2.6), we split the boundary into d^\^ and \ Sfiin. This yields 



\_d_ 
2di 



B^ 

n 



dxdy ~^ ^ J (curl(B)^ dxdy 



< a / B^ dxdy - ;^ ( / + / ) B^ (u • n) ds 

^ Ionia Jdn\anin J 



2 



/ + / ) (b • (curl(B) X n)) ds. 

Jdn,„ Jdn\dnin I ^ ' 



Rearranging the above relation and applying mixed boundary conditions (2.6), 
remembering that /? < — l/(2e), we get 

\^ \ B'^dxdy + e [ (curl(B)) dxdy 
'2 ot Jq^ Jq\ / 

<a [ B^ dxdy - [ ( - + e/3 ) B^ (u • n) ds. 

Jn Ja^in V2 / 

From the above relation, after integrating in time then we have the required result. 

□ 



3. Semi-discrete Schemes 

To simplify the treatment of the boundary terms we let the computational do- 
main be the unit square. A justification for this will be provided at the end of 
this section. 

The SBP finite difference schemes for one-dimensional derivative approximations 
are as follows. Let [0, 1] be the domain discretized with xj — jAx, j ^ 0, . . . , N — 1. 
A scalar grid function is defined as w = {wq, ...wn-i). To approximate dxW we use a 
summation- by-parts operator = P~^Qx, where is a diagonal positive N x N 
matrix, defining an inner product 



{v, w)p^ = v'^PxW, 
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1/2 

such that the associated norm ||ui||p — {w,w)p is equivalent to the norm lluijl — 
(Ax w^)^/^. Furthermore, for to be a summation-by-parts operator we 
require that 

Qx + = Rn — Ln, 

where Rn and Ln are the N x N matrices: diag (0, . . . , 1) and diag (1, . . . , 0) re- 
spectively. Similarly, we can define a summation- by-parts operator Dy = Py^Qy 
approximating dy. Later we will also need the following Lemma, proven in [ ]. 

Lemma 3.1. Given any smooth function u{x,y), we denote its restriction to the 
grid as u and let w be a smooth grid function. Then 

(3.1) \\Dx(uow) -~ no D^w\\p^ < C ||5a;u||ioo([o,i]) II^^IIp^ 

where {u o v)j = UjVj. 

Next, we move on to the two-dimensional case and discretize the unit square 
[0,1]^ using NM uniformly distributed grid points {xi,yj) = {iAx,jAy) for i = 
0, . . . , iV - 1, and j = 0, . . . , M - 1, such that {N - l)Ax = (M - l)Ay = 1. We 
order a scalar grid function w{xi,yi) column vector 

W = (wo,0, Wo,l, . . . , Wo,(J\/-l), Wi,o, ...... • , W(jv_i)^(Af_i))^ . 

To obtain a compact notation for partial derivatives of a grid function, we use 
Kronecker products. The Kronecker product of an Ni x N2 matrix A and an 
Ml X M2 matrix B is defined as the TViMi x N2M2 matrix 

^ OiiB . . . oipf^B 

(3.2) A(^B= : : 

yo-NiiB ... ajViWaSy 

For appropriate matrices A, B, C and D, the Kronecker product obeys the following 
rules: 

(3.3) {A ® B){C ®D) = {AC ® BD), 

(3.4) {A®B) + {C®D)^{A + C)®{B + D), 

(3.5) {A®Bf ^{A^ ®B'^). 

Using Kronecker products, we can define 2-D difference operators. Let /„ denote 
the n X n identity matrix, and define 

i>X ^ Da; (E) Im, dy^lN<E)Dy. 

For a smooth function w{x,y), {'i>xw)i.j ~ dxw{xi,yj) and similarly {Dyw)ij « 
dyw{xi,yj). 

Set P = Px (E) Py, define {w, v)p = w'^Pv and the corresponding norm ||w||p — 
1/2 

{w,w)p . Also define 72. = RnE)Im, ^ — LnE)Im,1^ — InE)Rm and I? = InE)Lm- 
For a vector valued grid function V — {V^, V'^), we use the following notation 



and so on. In the same spirit, the P inner product of vector valued grid functions 
is defined by (V, W)p = {V\ W^)p + {V^,W^)p. 
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Remark 3.1. Note that the Kronecker products is just a tool to facilitate the no- 
tation. In the implementation of schemes using the operators in the Kronecker 
products we can think of these as operating in their own dimension, i.e., on a spe- 
cific index. Thus, to compute DxW, we can view w as a field with two indices, and 
the one- dimensional operator will operate on the first index since it appears in 
the first position in the Kronecker product. 

The usefulness of summation by parts operators comes from this lemma. 

Lemma 3.2. For any grid functions v and w, we have 

{v,dxw)p + {'OxV,w)p = v'^[{TZ-£){lN<»Py)]w 

(3.6) 

{v,'Oyw)p + {dyV,w)p =v [{U - T)){Px (g) Im)] w. 
Observe that this lemma is the discrete version of the equality 

V (d^w) dxdy + I j {dxv)wdxdy= v{l,y)w{l,y) - v{0,y)w{0,y) dy. 
J Jn Jo 

Proof. We calculate 

{v, V^w) = {P^ (g Py) {P~^Qx ® Im) w 

= V^Qx (gl PyW 

= -v'^QI PyW + v'^iQx + Ql) ® PyW 

= -{P-^Qx (g iMvf [Px (S)Py)w+ {Rn - Ln) O PyW 

= - i^xv f (Px (SPy)w + v'^ {n - C) {In O Py)w. 
The second equality is proved similarly. □ 

For a vector valued grid function V, we define the discrete analogues of the curl 
and the curl(curl) operators by 

mrl(V) =5^y2_j,^yi and 



cutf (V) 



-dyy d:,y \ fV^\ _ f-^yyV^+dxyV^ 



^xy -^xxj \V^J \^xyV^-dxxV^ 



where dxy = <}x^y = ^y<^x etc. Before we define our numerical schemes, we collect 
some useful results in a lemma. 

Lemma 3.3. 

(V,cure(V))p = ||mrl(V)||^ 

(3.7) + [V^ [{U - V) {Px Im)] cut[(V) 

- [VY [(^ - >C) {In Py)] curl (V) . 

If u is a grid function, then 

(V, u o j),v)p = ^v^[{n - c){In ® Py)] [u o V) 



(3.8) 



+ - (u o dxV -dxiuoV), V)p , 



(V,«o5^V)p = lY^m-V){Px®lM)] {uoY) 

+ l{uoT,yV--Oy («oV),V)p. 
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Proof. To prove (3.7), 

(V, CUrP (V))p = {V\-dyyV^ + Dy^V^)p + {V^,d,yV' - d^^V^)^ 

= - {dyV\ -dyV' + ^^V^)p - (0.^', ^yV - O.V^')p 

+ [(^ - (P-:^ ® hi)] i-^yV + 0,^2) 

+ {VY [(7^ - C) {In ^ Py)] {dyV^ ~ 0,1/2) 

= (cur[ (V) , curl (V))p + (v'f [{U - V) (P, ® Im)] curl (V) 

- {V''f[{n-C) [In ® Py)]tuz{{\) . 

To show (3.8), first note that since P is diagonal, (uo9V, V)p = (OV, uoV)p. We 
use Lemma 3.2 to calculate 

(uoc),V,V)p ^ (c),(uoV),V)p + (iioc),V-5, (woV),V)p 

= - (u o V, c):,V)p + (7^ - C) {In ® Py) [u o V) 

+ (u o Xs^M - Oj, (u o V) , V)p 
= - (u o O^V, V)p + (7^ - C) {In ® Py) {u o V) 

+ (uoc)xV-0, (uoV),V)p. 

This shows the first equation in (3.8), the second is proved similarly. □ 

Now we are in a position to state our scheme(s). For = 1 or 2 we will use 
the notation for both the grid function defined by the function u^{x,y) and for 
the function itself. Similarly, for the boundary values, we use the notation h and g 
for both discrete and continuously defined functions. Hopefully, it will be apparent 
from the context what we refer to. 

The differential equation (1.5) will be discretized in an obvious manner. We 
incorporate the boundary conditions by penalizing boundary values away from the 
desired ones with a 0{1/Ax) term. To this end set 

B = [{P-^ ® Im) {^cC + ^nn) + {In ® Py') {^vV + S^^W)] , 

where S]£, Sp and Yiu are diagonal matrices, with components {o'c)jj ordered 
in the same way as in ((3.2)) (and similarly for the other penalty matrices), to be 
specified later. Furthermore, the following form of the penalty paramaters will be 
convenient: 

(3.9) B = B' + eB" 

and similarly for I]£, ac etc. 

With this notation the scheme for the differential equation (2.1) with boundary 
values B(a;,y) = Q{x,y) reads 

(3.10) Vt + o i)^ V + o 0,y V - CV + ecurP ( V) = S( V - g) , t > 0, 
while V(0) is given. Here C denotes the matrix 

^ _ f-'^yU^ ^yu' \ 
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Theorem 3.1. Let V he as solution to (3.10) with g = 0. If the constants in B is 

chosen as 

(3.11) 

{<^1Z)N~1,] < ^ , {<7c)0,j < 2 ' (.'^u)i:M-l < ^ , 

and (crp)i,o < ^ , 

and all other entries are 0, then 

(3.13) l|V(Ollp<e^*|lV(0)||^, 

where u^'^ = (u'vO), m''^ ~ (u'AO), fori = 1,2 andc is a constant depending onu^, 
V? , and their derivative approximations, hut not on N or M . By construction of the 
SBP operators p = p^ = pff_^ = Pq = Pm-i, where = Axdiag {pi,-.- ,P%-i) 
and py = Aydiag {p^, . . . ,p^/_i). 

Proof. Set E{t) = (V, V)p. Taking the P inner product of (3.10) and V, we get 
1 d 
2~dt 

+ (V,CV)p + (V,6V)p 

Using Lemma 3.3 we get 
l|i. + .||cu..(V)||^ 



-E + e (V, cud2(V))p - - (V, o c),V)p - (V, o 0^ V)p 



e(yi) [{U-V){P^®lM)]cux{{W) + e{V^) [{U - C) {I n ® Py)] tuzi [Y) 
- (V, o l),V)p - (V, o c)^V)p + (V, CV)p + (V, 6V)p 
(yi)^ [(Z^ - V) (P, $5 /m)] curl (V) + e (V^^)^ [(7^ _ £) (j^ ^ p^)] ^ucl (V) 



-e 



" [(7^ - £)(/Ar ® Py)] (u' o V) - [(Z^ - 2?)(P, ® hi)] {u^ o V) 
- i (yi o O.V - 0,(7/1 o V), V)p - i (^2 ^ _ j,^(^2 o V), V)p 
+ (V,CV)p + (V,(6' + eS")V)p. 
Note that by (3.1), 

I (lii o O.V - 0,(^1 o V), V)p| < c ||V||^ , 

(3.14) |(m2o0j,V-£)j,(u2oV),V)p| <c||V||^, 

|(V,CV)p|<c||V||^, 

for some constant c depending on the first derivatives of and u^. Using the 
conditions (3.11) we arrive at 

^|i. + .||cu.((V)||^ 

<cE~e {V^ [{U - V) (P, ® Im)] curl (V) + e {V^f [(7^ - £) (/^ ® Py)] curl (V) 
+ e(V,S"V)p. 
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Next, for any grid function w (with components as in ((3.1))), we have 

Ml = Aypy J2 ^^vl 



Similarly 



j=0 1=0 

N-l 

i=0 



M-l 

|w||p > Axp ^ Ayp^j (wl 

3=0 



- w: 



(3.15) 



2 . Axp 



Combining this we find that 

2 'w^(^{lN^Py)C+{lN^Py)n)w 

-T ^w^(((Px ® Im)V + (P, Im)u)w. 
We also compute 



(3.16) 



(u), Bw) 



+ {Px ® Im)^vV + [P^ ® Im)^uU)w. 
Using (3.15) for w = cur[(V), and (3.16) for = V we find 



dt 



E <cE 



cud(Vf {{In «) Py)C + {In (E> Py)n)cud(V) 
+ ^cur[(V)^((P, ® /m)I? + {Px ® lM)U)cml{V) 

+ {V^f [{U - V) (P, ® Im)] cur[ (V) - [V^f [(7^ - C) {In <E> Py)] curl (V) 

+ V^(-^2(^JV ® Py)C - a'^{lN ® Py)n - a'i{P, ® Im)V - a'^{P, ® /m)W)V 
=: -eA 

Choose the remaining penalty parameters as in (3.12) and write 
A>Pf (cur((V), cur[(V))^ + (V^, cur[(V))^ - (V^ V^) 



Pdx 



C 

;cur[(V), cur[(V))^ - {V\ cur[(V))^ - a'^ {V\V^) 



+ ^ (cur((V), cur[(V)>^ - {V\cux\{Y))^ a'i {V\V')^ 

Pdy 

2 

1/1 
2 



(curl(V), cur[(V))^ + {v\cml{V% - a'l^ {V\V^) 



lu 



1 



^JVdx 



2pAx 



1 / 1 



2 \ ^JVdx 



1 



1 



IVdx 



/p^zux{{Y) 



n 



'n 



2pAxJ ^ ' 
1 \ , , 
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>o, 

where 

{v,w)^ = v^{In ® Py)Tlw and so on are positive semi-definite forms. 

Furthermore, we used the notation pdx = pAx and pdy = p^y- Summing up, we 
have shown that 

j^E{t) < cE{t), 

and the result follows by Gronwall's inequality. □ 

The scheme for the mixed boundary conditions reads 
VtW o O^V + u^o dyV -CV + ecuc[2(V) 

where h is the desired value of cur[(V) on the boundary. 

Theorem 3.2. IfV is a solution of (3.17) with h = g = 0, and B{= B') is chosen 
so that (3.11) holds, then 

(3.18) \\Y{t)\\l + e f ||cur[(V(s))||^ ds < e^* |iV(0)||^ , 

Jo 

where c is a constant depending on S^,, Oj, and the derivatives of and v?' , hut not 
on N or M . 

Proof. The proof of this theorem proceeds as the proof of Theorem 3.1. Note that 
we are subtracting the boundary terms coming from (V, curP(V))p, so that we 
do not need to split B. The other terms are estimated as before, and we get the 
inequality 

^i? + e II cur[(V) II ^ <ci?, 

which yields the stability result. □ 

Remark 3.2. We have assumed a constant resistivity co-efficient e in the above 
discussion. However, in many practical applications, the co-efficient of resistivity 
can vary in space. In such cases, our theoretical results hold provided that the 
resistivity co-efficient is uniformly bounded away from zero. 

The analysis has been carried out on a Cartesian equidistant grid on the unit 
square. However, this is not a restriction as problems on general domains may 
be addressed using coordinate transformations. The stability proofs will hold as 
long as the norm matrices {Px.Py, Pz) are diagonal. SEP finite-difference schemes 
with diagonal P matrix have a truncation error of 2p in the interior and p near the 
boundary resulting in a global order of accuracy /convergence rate of p+ 1. (See [9] 
for further details.) 

The resistive magnetic induction equations include diffusive terms. Those are 
discretized by applying the first-derivative operators twice. This results in a trunca- 
tion error of p — 1 near the boundary for the diffusive terms. However, thanks to the 
energy stability of the scheme, the global convergence rate and order of accuracy 
remains at p -f 1 (see [ ]). 
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4. Schemes in Three-dimensions 

In this section, we are going to write down the three-dimensional version of the 
finite difference scheme for the equation (1.4). To begin with, we discretize unit 
cube [0,1]"^ using uniformly distributed grid points {xi,yj,Zk) — (iAx, jAy,kAz) 
for z = 0, . . . ,iV-l, j = 0, . . . ,M-1, and fc = 0, . . .,K~l such that (A^-l)Aa; = 
(Al — l)Aj/ = (K — l)Az ~ 1. We order a scalar grid function w{xi, yi, Zk) — Wijk 
as a column vector 

W = (wo,o,0, WQ,Q,1, ■ ■ ■ , Wo^o,(K-l): ■"^0,1:0, Wo,l,l, . . . ,Wqxm-i),(k-i), ■ ■ ■ , 

Wlfl,0: U'(Ar_i),(M-l),(K-l) j • 

As before, let /„ denote the n x n identity matrix, and define 

dx ^ Dx <E) Im <E) Ik, dy ^ In ^ Dy (g) Ik , D^, ^ In ® hi ® D^. 
Set P — Px®Py®Pz, define {w,v)-p — w'^Pv and the corresponding norm ||w||p ~ 

1/2 

{w, w)-p . Also define TZ — Rn ® Im ® Ik, ^ — Ln ® Im ® Ik,U — In ® Rm ® Ik 
and V — In ® Lm <X) Ik, A — In ® Im ® Rk, B = In ® Im ® Lk- 

For a vector valued grid function V — (y^, , V^)^ we use the following notation 

\^xV^ 

and so on. In the same spirit, the P inner product of vector valued grid functions 
is defined by (V,W)p = [V^ ,W^)p + {V^ ,W'^)p + {V^,W^)f. Finally, we set 

S = [{Px^ ® Im ® Ik) {^c^ + ^n'R-) + {In ® Py^ ® Ik) (Spl? -I- ^ul^) 

+ {In ® Im ® P~^) {^aA + Ee^)], 

where E^, T,n, E-p,Eiy,I]_4 and Eg are diagonal matrices. With these notations 
above the scheme for the differential equation (f .4) with boundary values B(x, y, z) = 
g(a;, y, z) reads 

(4.1) Wt+u^o X)xY + u^o DyV + u^o O^V -CY + ecuc[2(V) = S{V - g), t> 0, 

while V(0) is given. Here C denotes the matrix 

'()yU^ — dzU^ dyU^ 5zU^ 

C = I dxu'^ —^xU^ — fzW^ fzU^ 



and 




curl" (V) = I Oz {^yV^ - OzV^' - dx {dxV^ - ^yV^ 



y ' 

2\ 



^xV'')~dy {i)yV^-l>zV 

Remark 4.1. Note that the stability result for the three-dimensional scheme can be 
achieved along the same way as in Theorem 3.1 by choosing the penalty parameters 
as 

(4.2) 

, ^ {^,yj,Zk) IIS , u^'^{'i,yj,Zk) , , (xj,i,Zfc) 

\^n)N-i,j,k S ^ , (,o-£joj,fe S ^ , (cru)i,M^i,k S ^ , 

, , , U^''^ (Xi,0, Zh) , I . V?'^ {Xi.yj.l) , , M'^'+(X;, w,-, 0) 
{^v)^fl.k < < 2 ' " ^ "^''^^^ 2 
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and all other entries are 0. 



5. Numerical Experiments 

The SBP-SAT schemes (3.10) have been tested on a suite of numerical experi- 
ments in order to demonstrate their effectiveness. We have used the second-order 
(first-order) accurate and the fourth-order (second-order) accurate SBP operators 
in the interior (boundary) . From the results of [ ] , these operators result in overall 
second and third order accurate discretizations of the equations. Henceforth, the 
second (first)-order accurate SBP scheme will be denoted as SBP2 and the fourth 
(second)-order accurate SBP scheme will be denoted as SBPA after their orders of 
accuracy in the interior. Time integration is performed by using a standard second- 
order accurate Runge-Kutta scheme. (Using a higher-order Runge Kutta scheme 
didn't affect the quality of the computational results.) 

Numerical experiment 1: To begin with, we consider a divergence free velocity 
field u(a;. y) — (— y, x)^ and a slightly modified form of (2.1) given by 



(5.1) Bt + AiB^ + AsBy - CB = e 



iiB'^)xx - {B^)xy) 



where the forcing function T is given by, 
(5.2) 

/i = 160e(y-0.5sin(t)) [-4 + 40{(a; - 0.5 cos(i))2 + (t/ - 0.5 sin(t))2}] e^^^'^ 
/2 = -160e(y - 0.5cos(i)) [-4 -I- 40{(x - 0.5cos(i))2 ^ {y - 0.5 sin(t))2}] e^^^*', 
with A: 

A(t) = -20{(a;cos(i) + ?/sin(t) - 0.5)^ -I- (-x sin(<) + y cos(t))^} 

We note that it is straightforward to extend the stability results of the previous 
section to SBP-SAT schemes for (5.1). The forcing term is evaluated in a standard 
manner. The forcing function in (5.1) enables us to calculate an exact (smooth) 
solution of the equation given by, 

(5.3) B(x, t) = i?(t)Bo(i?(-t)x), 

where i?(t) is a rotation matrix with angle t. Note that this exact solution rep- 
resents the rotation of the initial data about the origin. In fact, (5.3) is also an 
exact solution of (5.1) with both the forcing term T and the resistivity e set to 
zero. Hence, this example follows from a similar example for the inviscid magnetic 
induction equations considered in [11, •">]. 

For initial data, we choose the divergence free magnetic field: 

(5.4) Bo(x,j;)=4(^-yi)e-M(-V2)^+^^), 

and the computational domain [— 1, 1] x [— 1, 1]. In this case, the exact solution is a 
smooth hump (centered at (1/2,0) and decaying exponentially) rotating about the 
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origin and completing one rotation in time t — 2tt. The hump remains completely 
inside the domain during the course of the rotation. Since the exact solution is 
known in this case, we use this solution to specify the data for the boundary con- 
ditions (2.3) or (2.6). The above setup is simulated using the SBP2 and SBP4 
schemes. Using Dirichlet or mixed boundary conditions led to very similar re- 
sults. Hence, we present results only with the mixed boundary conditions (2.6) in 
this case. The time-integration was performed with a second-order Runge-Kutta 
method at a CFL number of 0.5. The resistivity e = 0.01 was used. Also we have 
used the value of the penalty parameters as mentioned in Theorem 3.1, for instance 
we used cc = i^'r + ec^ = 



2 



2pAx 



and so on. We plot the P norm of the 

magnetic field: B = y/ (B^)^ + (B^)^ , at times t — n (half rotation) and t — 27t 
(full rotation) for both the SBP2 and SBP4 schemes in figure 5.1. As shown in 



■ ■■ 

I 



I 



(a) SBP2, half rotation 



(b) SBP2, full rotation 



■■■■■I ■ ■■ ■■■■■ ■■ in ^ 



H 



■ I 



(c) SBP4, half rotation 



(d) SBP4, full rotation 



Figure 5.1. Numerical Experiment 1: ||B|| — + for the 

SBP2 and SBP4 schemes on a 100 x 100 mesh. 



the figure, both SBP2 and SBPA schemes resolve the solution quite well. There 
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are very few noticeable differences between the second and fourth order schemes 
at this resolution. The shape of the hump is maintained during the rotation. A 
quantitative view of the results is presented in Table 5.1 where 

(5.5) relative percentage error = — — j- x 100, 

B„ti„i is the numerical approximation and is the exact reference solution and 
— (Ax w^)^/^. The errors are computed at time t — 2tt (one rotation) on 



Grid size 


SBP2 


rate 


SBP4: 


rate 


40x40 


2.1e-l 




1.6e-2 




80x80 


5.7e-2 


1.9 


l.le-3 


3.9 


160x160 


1.3e-2 


2.1 


l.le-4 


3.3 


320x320 


3.1e-3 


2.0 


1.3e-5 


3.1 


640x640 


7.5e-4 


2.0 


1.6e-6 


3.0 



Table 5.1. Relative percentage errors in for ||B|| at time t = 
27r and rates of convergence for numerical experiment 1. 

a sequence of meshes with both the SBP2 and SBPA schemes. The results show 
that the errors are quite low, particularly for SBPA and the rate of convergence 
approaches 2 for SBP2 and 3 for SBPA. This is consistent with the theoretical 
order of accuracy for SEP operators (see [i'*]). The very low values of error with 
SBPA suggest that one should use high order schemes to resolve interesting solution 
features. 

Another feature of numerical solutions of equations (5.1) is the behavior of diver- 
gence of the magnetic field. Note that both the initial data and the forcing function 
(5.2) are divergence free. Hence, the divergence of the exact solution of (5.1) should 
remain zero for all time. However, as remarked before, we don't attempt to preserve 
any particular discrete form of divergence. Hence, numerical divergence errors are 
an indicator of the performance of the schemes. We define the discrete divergence 
operator: 

divp(y) = 'o^v^ + dyV^. 

This corresponds to the standard centered discrete divergence operator at the cor- 
responding orders of accuracy. The divergence errors in P and rates of convergence 
at time t = 2tt for the SBP2 and SBPA schemes on a sequence of meshes are 
presented in Table 5.2. From Table 5.2, we conclude that although the initial di- 



Grid size 


SBP2 


rate 


SBPA 


rate 


20x20 


8.9e-l 




4.5e-l 




40x40 


3.9e-l 


1.1 


8.0e-2 


2.9 


80x80 


l.Oe-1 


2.0 


4.2e-3 


3.8 


160x160 


2.7e-2 


1.9 


5.0e-4 


3.0 


320x320 


9.5e-3 


1.5 


8.0e-5 


2.6 



Table 5.2. Numerical Experiment 1: Divergence errors in P and 
rates of convergence at time t — 2n. 

vergence is zero, the discrete divergence computed with both the SBP2 and SBPA 
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schemes is not zero. However, the divergence errors are very low in magnitude even 
on fairly coarse meshes and converge to zero at a rate of 1.5 and 2.5 for SBP2 and 
SBPA scheme respectively. A simple truncation error analysis suggests that these 
rates for the SBP2 and SBP4 schemes are optimal. Finally, we emphasize again 
that the quality of the solutions are good and the convergence rates did not suffer, 
despite the scheme not preserving any form of discrete divergence. 

Numerical Experiment 2. In the previous numerical experiment, the hump (rep- 
resenting the interesting parts of the solution) was confined to the interior of the 
domain. A more challenging test of the boundary closures is provided if the hump 
interacts with the boundary. We proceed to test this situation by considering (5.1) 
in a domain [0, 1] x [0, 1] with exactly the same initial data, resistivity and forcing 
function as in the previous numerical experiment. The exact solution (5.3), being a 
rotation about the origin, now exits the domain first at the lower boundary (includ- 
ing a corner) and enters the domain through another part of the boundary during 
the course of a single rotation. We study this interaction by simulating (5.1) with 
the SBP2 and SBP4 schemes. We consider (5.1) with Dirichlet boundary condi- 
tions. The boundary data are calculated by evaluating the exact solution (5.3) at 
the boundary. 

The norm of B after half a rotation and one full rotation is plotted in figure 
5.2. The figure shows that both the SBP2 and SBP4 schemes resolve the solution 
quite well and maintain the shape of the hump. Furthermore, the boundary inter- 
actions are resolved in a stable and accurate manner, showing that the choice of 
boundary closures was proper. The errors in are shown in tables 5.3 and 5.4 and 
we observe that the errors are quite low and the correct rates of convergence are 
obtained. The results were very similar to those obtained in numerical experiment 
1. 



Grid size 


SBP2 


rate 


SBP4 


rate 


20x20 


5.5e-2 




1.5e-2 




40x40 


l.Oe-2 


2.4 


1.9e-3 


2.9 


80x80 


2.3e-3 


2.2 


1.6e-4 


3.5 


160x160 


5.4e-4 


2.0 


1.9e-5 


3.1 


320x320 


1.3e-4 


2.0 


2.5e-6 


3.0 



Table 5.3. Relative percentage errors in and rates of conver- 
gence for numerical experiment 2. 



Grid size 


SBP2 


rate 


SBP4 


rate 


20x20 


5.6e-2 




9.8e-2 




40x40 


3.3e-2 


0.8 


3.7e-2 


1.4 


80x80 


9.3e-3 


1.9 


2.3e-3 


4.0 


160x160 


2.1e-3 


2.1 


2.5e-4 


3.2 


320x320 


7.4e-4 


1.5 


4.1e-5 


2.6 



Table 5.4. Divergence errors in P and rates of convergence at 
time 1 = 2^ for numerical experiment 2. 




Figure 5.2. Numerical Experiment 2: |B| = ^Bf + for the 
SBP2 and SBP4 schemes on a 100 x 100 mesh. 

This shows that the SAT technique for imposing boundary conditions weakly 
works very well even with complicated boundary data. 

Numerical Experiment 3. In the first two examples, we verified the accuracy 
of our schemes by manufactured solutions, using specific forcing functions. Next, 
we compute solutions of the magnetic induction equations in its original form (2.1) 
(without any forcing) to illustrate the role of the resistivity, e, in driving the dy- 
namics. We use the same velocity field and initial data as in the previous two 
numerical experiments. The domain is [—1,1] x [—1,1]. We compute solutions 
till t = 2tt with two different values of resistivity. The results with e — 0.05 and 
e = 0.001 are shown in figure 5.3. In the absence of exact solutions for B, we can 
only compare qualitative features in this case. For low resistivities like e = 0.001, 
the problem is very close to its inviscid version and the hump doesn't show much 
distortion during the rotation. The results in this case are very similar to the ones 
for the inviscid magnetic induction equations presented in [ ] . The SBP4 scheme is 
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(c) SBP4, e = 0.05 (d) SBP4, e = 0.001 

Figure 5.3. Numerical Experiment 3: |B| = ^Bf + for the 
SBP2 and SBP4 schemes on a 100 x 100 mesh. 



shghtly sharper (and hence more accurate) than the SBP2 scheme. Taking a higher 
value of resistivity e — 0.05, the viscous term starts playing an important role in 
the dynamics and the hump is expected to be smeared out. This is clearly shown 
in figure 5.3. The results of SBP2 and SBP4 schemes are very similar in this case. 

Since we do not enforce the divergence constraint excactly, we can use divergence 
errors as a quantitative measure. We know that the initial div(B) = and it should 
remain so during the computation. In tables 5.5 and 5.6, we display the divergence 
errors and their convergence rates for two different values of e. The convergence 
rates are as expected and we note that the errors are much lower for the higher 
order scheme. Thus, this experiment illustrates that both schemes are quite robust 
and efficient for different values of the resistivity. 
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Grid size 


SBP2 


rate 


SBP4 


rate 


20x20 


l.OeO 




7.3e-l 




40x40 


8.0e-l 


0.4 


1.2e-l 


2.6 


80x80 


2.7e-l 


1.6 


8.2e-3 


3.8 


160x160 


7.0e-2 


2.0 


l.Oe-3 


3.0 


320x320 


2.5e-2 


1.5 


1.7e-4 


2.6 



Table 5.5. Divergence errors in and rates of convergence at 
time t — 2t: for numerical experiment 3 with e = 0.001. 



Grid size 


SBP2 


rate 


SBP4: 


rate 


20x20 


7.3e-l 




1.4e-2 




40x40 


5.0e-2 


0.8 


2.2e-3 


2.7 


80x80 


l.le-2 


2.2 


1.7e-4 


3.7 


160x160 


2.9e-3 


1.9 


2.1e-5 


3.1 


320x320 


9.7e-4 


1.6 


3.4e-6 


2.6 



Table 5.6. Divergence errors in and rates of convergence at 
time t = 2t: for numerical experiment 3 with e = 0.05. 



6. Conclusion 

We have presented finite difference schemes for the magnetic induction equa- 
tions with resistivity. These equations arise as a sub model in the resistive MHD 
equations of plasma physics. We have shown that the symmetric form (1.4) of the 
magnetic induction equations with resistivity is well-posed with general initial data 
and both Dirichlet boundary conditions as well as mixed boundary conditions. 
SBP-SAT based finite difference schemes were designed for the initial-boundary- 
value problem corresponding to the magnetic induction equations with resistivity. 
These schemes were based on the form (1.4) and use SBP finite difference oper- 
ators to approximate spatial derivatives and an SAT technique for implementing 
boundary conditions. The resulting schemes are high-order accurate and shown to 
be energy stable. 

The schemes were tested on numerical experiments illustrating both the stability 
as well as high-order of accuracy. We also the use of the divergence errors as 
a measure of the accuracy of the solution. The results show that the SBP-SAT 
approach is a promising method to simulate initial boundary value problems for 
more complicated equations like the resistive MHD equations. 



Appendix 

For the sake of completeness, here we will present all the SBP operators used in 
the analysis. We consider second and third order accurate finite difference approx- 
imations. 
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First order accuracy at the boundary: The discrete norm P and the discrete 
second order accurate SBP operator P~^Q approximating ^ are given by 



P = h 



(\ 



\ 



\ 



P-'Q 



/-I 
1 

"2 



V 





-1 



1 

2 



We have used the operator P ^QP to approximate 4 



dx2 



Second order accuracy at the boundary: The discrete norm P is defined as 



/ 48 



P = h 



\ 



59 
48 



43 
48 



49 
48 



V ■■ / 

The discrete difference SBP operator P^^Q approximating is given by 



P-^Q 



f 



59 
34 

_ 59 
86 






_ J_ 

2 



_59 
-_£8 

12 



__3_ 
34 


59 
86 



_ 2 
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